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Simulations of charge transport in amorphous semiconductors are often performed in microscopi¬ 
cally sized systems. As a result, charge carrier mobilities become system-size dependent. We propose 
a simple method for extrapolating a macroscopic, nondispersive mobility from the system-size de¬ 
pendence of a microscopic one. The method is validated against a temperature-based extrapolation 
[Phys. Rev. B 82, 193202 (2010)]. In addition, we provide an analytic estimate of system sizes 
required to perform nondispersive charge transport simulations in systems with finite charge carrier 
density, derived from a truncated Gaussian distribution. This estimate is not limited to lattice 
models or specihc rate expressions. 


I. INTRODUCTION 


Charge carrier mobility is the key characteristic of 
organic semiconductors. Experimentally, it can be ex¬ 
tracted from time-of-flight measurements^’^, current- 
voltage characteristics in a diode^’^ or field effect transis¬ 
tor®’®, pulse-radiolysis time-resolved microwave conduc¬ 
tivity measurements^, or other techniques®^^^. 

In amorphous organic materials the energetic land¬ 
scape sampled by a charge carrier can be rather rough, 
with the width of the density of states as large as 0.2 eV. 
As a result, charge transport in thin films becomes dis¬ 
persive, i.e., the extracted mobility varies with the film 
thickness^®^^®. Consequently, the intrinsic value of mo¬ 
bility is difficult to measure: for example, the film thick¬ 
ness has to be large enough in time-of-fiight experiments, 
imposing stringent requirements on the accuracy of mea¬ 
surements of transient currents. 

A similar situation is encountered in computer sim¬ 
ulations of charge transport in organic semiconductors. 
Here, both lattice and off-lattice models employ system 
sizes which are usually much smaller than those used in 
experimental setups. This leads to an artificial increase 
of the average charge carrier energy and, as a result, to 
overestimated values of the charge mobility^^^^®. 

To overcome the limitations imposed by small system 
sizes, a method based on a temperature-extrapolation 
procedure has recently been proposed^^. Its main idea 
is to simulate charge transport at a range of elevated 
temperatures. At high temperatures transport becomes 
nondispersive and one can then extract the nondisper¬ 
sive mobility at relevant (lower) temperatures from the 
mobility-temperature dependence, /i(T). This method 
relies on the analytical dependence of the mobility on 
the temperature derived for a one-dimensional system^^ 
with Gaussian distributed energies and Marcus rates for 
charge transfer (see Methods section), given by 
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( 1 ) 


In three dimensions Rq, a, and b are treated as fitting pa¬ 
rameters instead of the parameters derived analytically 


in one dimension. Hence, for three-dimensional trans¬ 
port, Equation (1) has to be validated for every particular 
system. It would therefore be useful to have an alterna¬ 
tive approach, which does not rely on the ad-hoc /r(T) 
function. This is the first target of our paper: solving 
the one-dimensional stochastic transport we derive the 
system-size scaling of mobility and benchmark it against 
the temperature-based extrapolation. 

Due to the filling of the energy levels in the tail of the 
density of states the charge density is known to have a 
strong impact on the mobility^®^®^. With increasing den¬ 
sity the energy per carrier is increased, leading to higher 
mobilities. On the other hand, in small systems mobili¬ 
ties are artificially increased. One can therefore presume 
that the error induced by finite-size effects will decrease 
with the charge carrier density. Hence, the second task 
of this manuscript is to provide a criterion for the system 
size required for nondispersive transport simulations at 
finite charge concentration. 


II. METHODS 

To perform mobility simulations, we use the Gaussian 
disorder modeP®’®^^®^, i.e., we assume that molecular 
sites are arranged on a cubic lattice and that site ener¬ 
gies, €i, follow a Gaussian distribution, /(e) = l/cr^/^x 
exp (—e^/2cr^), where tr denotes the energetic disorder. 
We assume a mean of zero throughout. We use the Mar¬ 
cus expression for charge transfer rates®®^®^ 

(Acij -|- qF ■ Vij -|- Xij) 
^XijkBT 

(2) 

for transitions j ^ i, where Ae^ = e^ — ej is the site 
energy difference, q is the charge, F is an external field, 
Tij = Ti — Tj is the distance between two sites, Xij is 
the reorganization energy and Jij denotes the electronic 
coupling. As a simplification, we assume a constant reor¬ 
ganization energy, Ay = A, a constant transfer integral, 
Jij = J, and a lattice spacing of |ry| = a. Further, T 
is the temperature and /cb the Boltzmann constant. The 
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Marcus rates allow to link the mobility to the chemical 
composition of organic semiconductors^^’^®^"^^. 

Charge-charge interactions are modeled by an exclu¬ 
sion principle"^®, i.e., each site can be occupied by only 
one charge carrier at a time. As a result, the equilibrium 
site occupation is given by Fermi-Dirac statistics^'^ 


P(e) 



( 3 ) 


where the Fermi energy, ep, is implicitly determined by 
the number of charges in the system through 


p{e)f{e)de = n. 


( 4 ) 


Here n is the charge carrier density, i.e., the number of 
charges divided by the number of sites. The average en¬ 
ergy per charge carrier, Cc, is then given by 


_ /.!lep(e)/(e)de 
/^P(e)/(e)de ’ 


Note that in the limit of zero charge carrier density or 
for high temperatures the Fermi-Dirac distribution can 
be approximated by the Boltzmann distribution, PB(e) = 
exp (—e/ZcbT), which yields Cc = —a^/k^T. 


III. SCALING RELATION 

A. Derivation 

We now derive and test the system-size dependence of 
the charge carrier mobility, /i(iV), in the limit of zero 
charge carrier density. The derivation is based on the 
model of a one-dimensional chain of length N with Gaus¬ 
sian distributed, uncorrelated energies, and hopping tak¬ 
ing place only between adjacent sites according to the 
Marcus rates. Equation (2). An electric field of strength 
F = |F| is applied in the direction of the chain. We will 
require the mean velocity 


vn{F) = {NK (6) 

where (•) denotes the average over the energetic disorder, 
i.e., {g (e)) = / def (e) g (e). Here we have approximated 
the mean rate by the inverse mean first passage time 
(tat), which can be calculated more readily. 

For completeness, we now present a detailed derivation 
of (tat). At steady state conditions for a given realization 
of the disorder, the mean first passage time to traverse 
the chain starting at * = 1 reads^® 


tn 


N-l i 


EE 


1 


^k,k+l 

Wfc+l,fc 


The rates fulfill detailed balance, /iop 
exp [—j3 (ci — Ck — f {i — k))] , which leads to 


n 

j=k 


^3,0 + -^ _ „-/3/(i-fc) 




exp<| /3^(ej+i -ej) 
j=k 


= g-/3/(*-fc)-|-/3(ei-efe) 


( 8 ) 


where / = qFa and /3 = Xjk^T. After shifting i-kt-d k, 
this can be rewritten to 


Tn 


N-l ^ i-1 

g-/3/fc+/3(ei-£i_fc) 

i=l 



( 9 ) 


For Gaussian distributed energies the exponential, e'^*, is 
log-normal distributed with mean F’ . Since, further¬ 
more, energies are uncorrelated, (eiCj) = cr^Sij, for fc > 0 
we have 

^g-/3Q-fc) (10) 

leading to 


(tn) = E ( 
2=1 ' 


N-l 


= E 


2=1 


^2+1,2 
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l + g5(/3'^)'^g-/3/fe 
/l-l 
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(geometric series with 2 = < 1). We can split the 

average because the first term only involves i and i -I- 1 
and the second term sites < i. The first term becomes 

/(/) = /j—\ = — Ldu+P/{i\)(u+^-M+x'f\ ^ (^ 2 ) 
/ Wo \ / 


where wq = 2t: AnXk^T is the prefactor of the rates 

Equation (2) and A' = A — / is the shifted reorganization 
energy. It is convenient to transform site energies as 

ei = e-^'5, €i+i=e + ^S. (13) 

The brackets above become 

+ (14) 


since integrals over other sites reduce to unity and the 
variable change has unity Jacobian. The integral over e 
can be evaluated straightforwardly, leading to 




+ |;(J + A')q, 



( 7 ) 


(15) 
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\/N temperature T (K) 

FIG. 1. (a) System-size extrapolation for energetic disorder of cr = 0.1 eV, external field ol F = 10® V/m, lattice spacing of 1 nm, 
transfer integral of J = 10~® eV and reorganization energy of A = 0.3 eV. (b) Validation using the temperature extrapolation. 
Large symbols denote the extrapolated mobilities, ^oo, summarized in Table I. 


Evaluating the second integral, we have 

+ 1/JA (l - {)'} (16) 

with 

A’ A-/3(t2 ^ 

Note that this result implies, in principle, an upper bound 
(/3 (t)^ < /3A for the disorder, cr, beyond which the mean 
waiting time diverges and, consequently, the approxima¬ 
tion in Equation (6) breaks down. The mean waiting 
time finally reads 

N-l 

(tn) 

i=l 

which constitutes our first central result. 
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(18) 


B. Mobility 


N-l 


thus obtain 

(T)v)=/( 0)^ [l + e^(/3<rf(*-l) 

= /(0)(1V - 1) [l + - 1) 

using 


z — z 

lim- = i — 1. 

z->l 1 — z 


For large N, 




(19) 


( 20 ) 


( 21 ) 


and the velocity decays as vn{0) ~ 1/N. One the other 
hand, evaluating the sum first, we obtain 

(™) = /(/) [(JV - 1) + 

( 22 ) 

which reduces to the same result as Equation (21) in the 
limit z —>■ 1 applying L’Hospital’s rule twice. 

For a finite field, by putting Equation (22) into Equa¬ 
tion (6) and taking the limit N ^ oo we get 


1 


— = /(/) 1 + 62 




1 - Z 




1 - z 


■ (23) 


With a few simplifications for sufficiently large N, we can 
thus write 


1 


1 


1 - 


1 1 


We first consider the limit z —> 1 corresponding to a 
vanishing field, F —>■ 0, before evaluating the sum. We 


'^N '^CO 
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T extrapolation 

N extrapolation 

200 K 

300 K 

600 K 

3.2 X 10”^® 

4.9 X lO”!'’ 

3.5 X 10"® 

2.7 X 10"^® 

2.8 X 10"“ 

3.3 X 10"® 


TABLE I. Extrapolated mobilities for the thermodynamic 
limit, fioo, from the temperature extrapolation method and 
the system size extrapolation method, all in m^V“^s“^. 

To leading order this results in 

fiv ~j > t^oo- (25) 


Specifically, for the zero-field mobility we obtain 



which is the same result that Seki and Tachiya have ob¬ 
tain in their Equation (6.7), cf. Equation (1). Going 
beyond their result and including the leading order cor¬ 
rection, we find for the mobility of the one-dimensional 
chain of N sites 

Mat = ficxj ^1 -|- (27) 

with c = |. Eor large N this result can be further sim¬ 
plified to In /tn ~ In Moo + ■ 

C. Numerical Test 



box size N 


FIG. 2. Convergence of the carrier energy with system size, 
N, for different carrier densities, n. Symbols are the results of 
random number experiments. Equation (28), while solid lines 
are the prediction of our analytic model. Equation (35). The 
dashed lines are the exact values for tc as N —>■ oo obtained 
from Equation (5). 


Similar to the temperature-based extrapolation, we 
now assume that Equation (27) also holds in three di¬ 
mensions, but with a different constant c. To test this 
assumption and to extrapolate the mobility value, we 
performed kinetic Monte Carlo simulations in cubic lat¬ 
tices from 8 X 8 X 8 to 50 X 50 X 50 sites, with Gaus¬ 
sian distributed energies and Marcus rates as described 
in the Methods section. The master equation for oc¬ 
cupation probabilities is solved using the variable step 
size algorithm^^d6,47_ charge mobility, ^ = d/rF, 
is evaluated using the charge trajectory, where d is the 
distance traveled by the charge along the field F during 
time r. The results are shown in Figure 1(a). One can 
see that the mobility (and its logarithm) scale indeed lin¬ 
early with the inverse system size, which is given by the 
number of sites in the system. 

The extrapolated mobilities, Moo = m(-^ oo): s-bo 
agree well with the temperature-based extrapolation, 
which is shown in Figure 1(b). Both methods are com¬ 
pared in more detail in Table I. 


IV. FINITE CARRIER DENSITY 

We now turn to the estimation of the error introduced 
by finite-size effects in systems with finite charge carrier 


density. A generalization of the approach of Sec. Ill to 
multiple interacting carriers is not straightforward since 
we are now faced with an exclusion process in the pres¬ 
ence of disorder, for which an analytical result for the 
mean first passage time is not available. Instead, we use 
the average energy per carrier, ec{N), as a figure of merit. 
This also makes the error estimate independent of the 
rate expression and the positional order of sites. Hence, 
it is also applicable to realistic morphologies and models 
with different charge transfer rates. 

We first performe a direct evaluation of €c{N) by draw¬ 
ing N random energies from a Gaussian distribution 
of width (T and calculating the ensemble average as 


ec(iV) 




(28) 


Results for a = 0.1 eV and T = 300 K are shown in 
Figure 2 (symbols) for different charge densities. This 
demonstrates that in hnite systems there is a significant 
deviation from Cc, especially at low charge carrier densi¬ 
ties. 

In order to obtain a closed-form expression for the 
finite-size error, we first note that the probability to draw 
an energy smaller than cq from a Gaussian distribution 
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Sec/ 

ec| < 5% 

Sec/ 

ec| < 0.1% 

\ fT(eV) 
n \ 

0.001 

0.01 

0.1 
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0.01 
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10® 
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10"® 

10® 

10® 

10® 

10® 

10® 

10® 

lO"'^ 

10® 

10® 

10® 

10® 

10® 

10" 

10"® 

10® 

10® 

10'‘ 

10® 

10® 

10® 

10“^ 

10® 

10® 

10® 

10® 

10® 

10® 

lO"! 

10® 

10® 

10® 

10® 

10® 

10® 


TABLE II. Necessary system size (order of magnitude) to en¬ 
sure that the relative error on the energy per charge carrier, 
Sec/sc, is smaller than 5% and 0.1%, respectively, assuming 
a temperature of 300 K. Errors are calculated using the dif¬ 
ference of the analytic estimate. Equation (35), to the exact 
value in an infinite system. Equations (5). 


f{x) reads 


T’ (e < eo) = f f{x)dx = F(eo), (29) 

J — CXi 

where F(x) is the cumulative distribution function 

FW = i + ierf(^). (30) 

The probability to draw an energy larger than eo is 
then given hy V {e > eo) = 1 — F{eo)- If we draw N 
independent energies, the probability that none of them 
will be smaller than eo reads 

IP (e, > eo , ^ = 1... A^) = (1 - F(eo))'^ . (31) 

The probability to find one value e < eo is then given by 

V {3i : Ci < eo) = 1 - V {ei > eo ,i = I... N) 

= l-(l-F(eo))^ (32) 

which is the cumulative distribution function for ep. The 
respective probability distribution for the minimum sam¬ 
pled value (MSV) function is obtained by differentiation 
and reads 


fusvix) = -N {F {x))^-^ fix). (33) 


The expectation value for the maximum sampled energy 
is given by Cmax = ~einin- With this model distribution 
function we obtain an estimate for the size-dependent 
average energy 


edN) 


tmax 

^min 


^max 

^min 


ep(e)/(e)de 
p(e)/(e)de ’ 


(35) 


which constitutes our second central result. This esti¬ 
mate is also shown in Figure 2 (solid lines) and is in 
good agreement with the values simulated directly. 

Hence, given the error, Sedn^a^N) = 
\ed'n,cr, N) — edn,a, N ^ oo)\, we can estimate the 
necessary system size N for our simulations. Such 
estimates are shown in Table II. As we have anticipated, 
large energetic disorder requires large system sizes, while 
for large charge densities one can use smaller systems. 

As of today, atomistically-resolved simulations can 
handle systems of approximately 5000 molecules. With 
coarse-grained models one can increase this number to 
about 10® molecules^®. In lattice models system sizes 
are up to 3 x 10® sites^®. Comparing this to Table II 
shows that for low carrier densities and for large values 
of energetic disorder, the necessary system size is compu¬ 
tationally still inaccessible, be it microscopic, stochastic 
or lattice models, and the extrapolation schemes from 
the previous section have to be used. For a sufficiently 
large charge carrier density or small energetic disorder, 
however, the error is within an acceptable range even for 
simulations in smaller systems. 


V. CONCLUSIONS 

To conclude, we have derived a system-size dependence 
of charge carrier mobility and provided a simple way of 
correcting for finite-size effects in computer simulations 
of charge transport in disordered organic semiconductors. 
We have also estimated the system sizes required for sim¬ 
ulating charge transport of carriers with a given error on 
the mean energy. Our results are general and are ap¬ 
plicable to different rate expressions as well as off-lattice 
morphologies. 
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